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The author reports on new high-fidelity simulations of charge carriers in the high-Tc cuprate 
materials using quantum Monte Carlo techniques applied to the first principles Hamiltonian. With 
this high accuracy technique, the doped ground state is found to be a spin polaron, in which charge 
is localized through a strong interaction with the spin. This spin polaron has calculated properties 
largely similar to the phenomenology of the cuprates, and may be the object which forms the Fermi 
surface and charge inhomogeneity in these materials. The spin polaron has some unique features 
that should be visible in X-ray, EELS, and neutron experiments. The results contained in this paper 
comprise an accurate first principles derived paradigm from which to study superconductivity in the 
cuprates. 


I. INTRODUCTION 

Understanding high temperature superconductivity in 
the cuprates is a long-standing and major challenge in 
condensed matter physics. The holes are the quasipar¬ 
ticles from which theories of the superconductivity are 
made, and in fact there have been many attempts to 
study doped holes using various computational and the¬ 
oretical techniques. The key question is the nature of the 
holes upon doping from the antiferromagnetic insulating 
state. The combination of calculations and experiments 
have resulted in a substantial amount of understanding 
of the holes in cuprates. An important concept is sponta¬ 
neous electron localization, proposed early on by authors 
such as Zhang and Ric^ and Emery and Reiteil^. These 
early proposals have been followed up by other authors 
suggesting mechanisms for electron localizatiorP^^. The 
literature on this su^ect is extensive and summarized in 
a number of reviews^^^. 

On the experimental side, a ngle resolved pho¬ 
toemission spectroscopy (ARPESj^®^ and quantum 
oscillatiorP^ measurements have elucidated the Fermi 
surface evolution as a function of doping, and gener¬ 
ally agree on a Fermi surface that starts near (7r/2,7r/2) 
in the crystallographic Brillouin zone and expands to a 
large pocket with doping. Optical spectroscop}^^ has 
found that, in addition to a Drude peak as the system be¬ 
comes metallic, additional absorption in the infrared part 
of the spectrum appears at around 1 eV. Scanning tun¬ 
neling microscopy (STM) has established that the s yst em 
is inherently inhomogeneouJI^®^. Finally, neutro iP, X- 
ra}EII, and STM technique^^ have found that these holes 
can arrange in stripe-like patterns. 

A major theoretical challenge is that traditional elec¬ 
tronic structure methods like density functional theory 
(DFT) suffer from severe errors in treatment of cor¬ 
relation. There have been studies using LDA+lP^^^ 
and hybrid functional^, and using quantum chem¬ 
istry techniques on cluster representation^^. Recently, 
DFT+DMFT has been applied as welP^. The pictures 
emerging from these calculations have many similarities 
and many differences as well. For example, while most 


techniques do correctly obtain holes occupying mostly 
the oxygen states, they disagree on whether the inter¬ 
action with the copper spins is ferromagnetic or anti¬ 
ferromagnetic in nature. There is a missing element in 
the theoretical techniques: a variational, explicitly corre¬ 
lated, ab-initio calculation of the bulk material, in order 
to separate proposals for hole states. 

In this article, I use highly accurate quantum Monte 
Carlo calculations of the ab-initio electronic structure 
of holes in the cuprates. This technique treats impor¬ 
tant short-range electron correlations accurately, which 
allows us to obtain a perspective on the effective low- 
energy electronic structure. It is also variational, which 
allows one to test different proposals for the hole state on 
equal footing. From these calculations, it appears that 
the most likely model for holes in cuprates is similar to 
one proposed by Emery and Reiteil^. I will show that this 
hole state is consistent with the measured Fermi surface 
and optical spectrum, and can accomodate the formation 
of stripes. 

The main purpose of this work is to establish the low- 
energy electronic structure of a hole in the cuprate mate¬ 
rial. At the level of accuracy currently feasible, quantum 
Monte Carlo techniques have the advantage that they 
are completely first principles and treat localized correla¬ 
tion very well; however, highly multiconfigurational wave 
functions are computationally out of reach given current 
capabilities. This work is thus concentrated on under¬ 
standing the ’building blocks’ of the low-energy physics 
of the cuprates and should serve to inform effective mod¬ 
els of their behavior. That is, if the ground state is highly 
multiconfigurational, then these single determinant-like 
states should be an important part of the correlated 
ground state. 


II. METHOD 

Fixed node diffusion Monte Carlo (FN-DMC) is a 
state-of-the-art method to calculate the electronic struc¬ 
ture of materials from first principles, and has recently 
been found to have high accuracy for the undoped 
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cuprate^^^^. Starting with a trial function I^t), the 
ground state is projected out by applying the imaginary 
time operator exp(—/fr). Exact projection suffers from 
the sign problem, which causes the method to scale ex¬ 
ponentially in the system size. The sign problem can 
be avoided by making the fixed node approximation, in 
which the zeros of the solution are constrained to be the 
same as the zeros of the trial wavefunction. This intro¬ 
duces a dependence on and the energy obtained is a 
variational upper bound to the true ground state energy. 
In this work, many different T^’s are considered to es¬ 
timate the best approximation to the ground state of a 
hole in the cuprates. Details of the calculations are very 
similar to ReP^ and are recorded in the Supplemental 
Information. 

The first principles Hamiltonian was used: 

H = + ( 1 ) 

i i<j ioL 

are electron indices, a is a nuclear index, and the 
nuclear-nuclear interaction has been omitted f or bre vity. 
Effective core potentials from Burkatski et al )^^ * ^^ * were 
used to eliminate the core electrons and give the form 
for Va- Density functional theory calculations were per¬ 
formed using CRYSTAlP^ to produce a starting Slater 
determinant, which was allowed to break spin symme¬ 
try to form localized moments. The determinant was 
varied by changing the starting magnetic order, which 
resulted in determinants with different arrangements of 
local moments (see Eigurej^, and by using hybrid den¬ 
sity functional theory calculations with a varying mixing 
parameter. Doped systems were simulated by removing 
one electron from a unit cell and compensating with a 
uniform background charge; the minimum cell size for a 
given doping percentage was used. The QWallP^ package 
was used to perform the quantum Monte Carlo calcula¬ 
tions. The Slater determinant was multiplied by a two- 
body Jastrow factor, which was variance optimized. Dif¬ 
fusion Monte Carlo was then performed using the Slater- 
Jastrow wavefunction as 

In diffusion Monte Carlo, there are a number of pa¬ 
rameters that determine the accuracy. All major pa¬ 
rameters have been checked to the highest degree pos¬ 
sible (Eigure[^. The timestep and finite size was var¬ 
ied, with no changes within stochastic errors. The nodes 
in the input Slater determinant were varied using dif¬ 
ferent hybrid density functional theories to generate the 
orbitals; this tuning adjusts between localized and delo¬ 
calized electronic structure. The minimum energy nodal 
structure was taken, which was always at 25% mixing. 
Einally, the dependence on the interlayer was checked by 
considering CaCu 02 and Ca 2 Cu 02 Cl 2 structures, with 
no change in the energy differences within stochastic un¬ 
certainties. As shown in Eigure[^ these parameters are 
converged. While the solution is not exact, in particular, 
long-range multiconfigurational character is not captured 
in this technique, these calculations are the highest ac¬ 
curacy ab-initio results for a hole in the cuprates. The 
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FIG. 1: Assessing errors in the calculations for a 2 x 2 x 1 
supercell, with x = 0.25. (a) The dependence of the rela¬ 
tive magnetic energies versus the density functional used to 
generate the nodes. All structures have a minimum within 
stochastic uncertainties at around 25% mixing, (b) At 25% 
mixing, the dependence on time step and on the interlayer. 
The physics is not qualitatively changed by these parame¬ 
ters. Stochastic uncertainties are approximately the size of 
the symbols. 


rest of the results in this paper will be for CaCu 02 with 
r = 0.02 Hartree”^, Slater determinants generated with 
the PBE functional at 25% exact exchange mixin^^^^^, 
twisted boundary conditions over real twists, and T- 
move^, unless otherwise indicated. 


III. RESULTS 
A. Total energy 

Eigure contains a summary of energetics for trial 
wave functions that differ in their magnetic ordering, for 
X = 0.25 and the 2\/2 x 2\/2 x 1 supercell. The ’clover- 
leaf’ shapes are the location of the Cu atoms, with oxygen 
atoms between them. The hole density is calculated by 
subtracting the doped charge density from the undoped 
AEM-ordered charge density. No matter the magnetic or¬ 
dering, the hole density is largely situated on the oxygen 
atoms, in agreement with X-ray experiments. Chang¬ 
ing the magnetic ordering affects the distribution of hole 
charge. 

An immediately striking result in Eigure is that the 
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FIG. 2: (a) FN-DMC energies of different magnetic orderings considered for a 2a/2 x 2a/2 x 1 unit cell of CaCu02. (b) 
Hole charge density obtained by subtracting the charge density of the x — 0.125 system from the x — 0.00 antiferromagnetic 
ordering, (c) Spin density of the corresponding orderings for x = 0.125. Both are calculated from the optimal single Slater 
determinant for that spin configuration for CaCu02 and projected onto the ab crystallagraphic plane. Correlations do not 
affect these pictures within statistical noise. In the density maps, blue is positive and red is negative. The density maps are 
normalized to the same value. 


flipped configuration, which has a single copper atom 
with spin reversed from the checkerboard AFM pattern, 
is lowest in energy. The flipped spin creates a region 
of five copper atoms with aligned spins, and a hole is 
attracted mostly to the oxygen atoms between the spin- 
aligned atoms. The ‘flip’ configuration is a spin polaron. 
The closely related compound, Ca2Ni02Cl2 does not ex¬ 
hibit this effect (Figure Q, so it appears to be a special 
feature of the cuprates. 


B. Excitation properties 

The gap was computed in FN-DMC by promoting an 
electron in the trial Slater determinant from the high¬ 
est occupied state to the lowest unoccupied state. In the 
case of Ca2Cu02Cl2, the PBEO ordering is incorrect, and 
the promotion was performed from the second-highest 
one-particle state. For the doped configurations, both 


spin channels were attempted and the lower energy ex¬ 
citation was chosen. In Figure the gap as a function 
of the twisted boundary condition in the supercell are 
presented for the lowest energy magnetic configuration 
at each doping level. Because of Brillouin zone folding, 
the twists available are limited to the ones shown. At 
X = 0.00, the minimal gap is 2.6(1) eV, a little above 
the experimental gap of about 2.0 eV for the undoped 
cuprates. The correction for calculating the gap at the 
F point is around 0.5 e\^, and so it is in good agree¬ 
ment with the experiment. Meanwhile, at x = 0.25, the 
gap clearly closes at O), which is in agreement with 

ARPES. 

At X = 0.125, it is possible that the gap closes near the 
F point (0,0,0), which, accounting for band folding for 
the 2^/2 X 2v^ x 1 unit cell, is where the ARPES Eermi 
surface is be located. Because of the larger supercell, the 
stochastic errors could not be reduced below around 0.1 
eV. Eor CaCu02, it appears that the gap is not quite 
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FIG. 3: Optical gaps in FN-DMC and PBEO calculated by 
promoting an electron in the Slater determinant, performed 
in the 2x2x1 cell for x — 0.00 and x — 0.25, and in the 2\/2 x 
2^f2 X 1 cell for X — 0.125. The twisted boundary condition 
of the wave function in the crystallographic primitive cell is 
labeled on the horizontal axis. 


closed; however, systems close to metal-insulator transi¬ 
tions often suffer from larger finite size effects. To check 
this, I also considered Ca2Cu02Cl2 at the F point, which 
has a larger c-axis direction. It appears that the gap is ei¬ 
ther zero or near zero at the F point for the Ca2Cu02Cl2 
model, which increases the likelihood that the gap is ac¬ 
tually closed at x = 0.125. 

There is also an excitation around 1 eV that ap¬ 
pears upon doping. This may correspond to new states 
seen in optical experiments at about that energ}^^. So 
the spin polaron configuration has excitation properties 
largely in agreement with those seen in experiment for 
the doped cuprates. Unfortunately the resolution is not 
high enough to comment on potential Fermi arcs. 


C. Polaron-phonon coupling 

While the focus of this article involves mainly the elec¬ 
tronic degrees of freedom only, since the lowest energy 
state involves localized charge density, one might expect 
coupling to the lattice. Indeed, in the hybrid DFT calcu¬ 
lations that also find a flipped ground state, the lattice 
reacts strongly to the presence of the hole, with relax¬ 
ations of approximately 0.1 A. In particular, the oxygen 
breathing mode is affected by the hole. This mode may 


be responsible for kinks in ARPES spectrsP^. These ef¬ 
fects warrant further investigation and will likely be im¬ 
portant for a full description of the spin polaron. 


D. Mechanism for the spin polaron 

The behavior of the hole seen in the FN-DMC results 
warrants some explanation. Actually, the unusual behav¬ 
ior of magnetism in the cuprates begins even before the 
material is doped. Even in the undoped regime, there is 
a ferromagnetic-like interaction between the oxygen and 
copper. This is in contrast to the normal picture, even 
confirmed recently in FN-DMC calculations on V02p^, 
in which the interaction between the ligand and transi¬ 
tion metal is antiferromagnetic-like. Once the material is 
doped, the copper spins change only by ^0.1 Bohr mag¬ 
netons, but the oxygen spin density is disrupted, as seen 
in Figure 

The unusual properties in the previous paragraph can 
be understood in terms of a pair of orbitals, diagrammed 
for a two copper unit cell in Figure with x = 0.00 and 
in the AFM ordering. The low-energy orbital is bonding¬ 
like between the oxygen and copper atoms, with more 
weight on one of the copper atoms due to the electron 
interaction. The second, higher energy, orbital is bond¬ 
ing for the spin minority copper atom and antibonding 
for the spin majority atom. When the spin density is 
calculated just from those orbitals, the low-energy or¬ 
bital exhibits the antiferromagnetic relationship between 
the copper and oxygen spins. However, the partial anti¬ 
bonding orbital has the opposite behavior and more spin 
density on the oxygen atoms. When the spin densities are 
summed, one obtains the ferromagnetic-like relationship 
between the oxygen and copper atoms. 

On doping, the partially antibonding orbital is depop¬ 
ulated in PBEO. Since this orbital is responsible for the 
coupling between adjacent spins, the coupling between 
copper spins changes from antiferromagnetic to ferromag¬ 
netic. Electron correlations are important in determin¬ 
ing that this orbital is depopulated upon introduction 
of a hole. In the FN-DMC calculations, we see that the 
spin up/spin down covariance on the copper atom around 
which the hole is centered decreases from -0.05 to -0.02. 
This decrease is due to a decrease in the double occu¬ 
pancy of that copper atom. In addition, the orbital is 
partially antibonding, which costs kinetic energy. The 
spin polaron, therefore, is stabilized by a balance between 
the kinetic energy and the interaction energy between 
electrons, which necessitates a correlated approach like 
FN-DMC. 


E. Charge and spin stripes 

PattersoiP noted that in hybrid DFT (B3LYP) calcu¬ 
lations, it was possible to form charge and spin density 
waves using spin polarons. In Figure a similar struc- 
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FIG. 4: Orbitals contributing to the spin density of a two 
formula unit cell of CaCu02. (a-c) Lower energy orbital that 
contributes to most of the spin density on the copper atoms, 
(d-f) Higher energy orbital that contributes to most of the 
spin density on the oxygen atoms. Blue is positive, red is 
negative. 


ture is presented, along with the Fourier transform of the 
charge and spin. Since the QMC calculations agree with 
DFT(PBEO) on the ground state density, it is likely good 
enough to analyze the properties of such a stripe system. 
The stripe structure in Figis in several ways quite close 
to that seen in neutron, x-ray, and STM experiments. It 
matches the Bragg peak in the charge density at (0.25,0), 
as well as a d-wave intracell density, as seen in STIVP^. 

If these stripes of spin polarons are the objects respon¬ 
sible for the stripes in the cuprates, then there should be a 
small peak in neutron diffraction at (1,0) and (0,1) due to 
the FM-like coupling between the copper atoms (Fig[^), 
in addition to AFM-like (1,1) peaks. This peak is a nec¬ 
essary prediction of this physics; if it is not present, then 
these objects cannot be responsible for the stripes. On 
the other hand, the Fourier transform of the charge den¬ 
sity should have small peaks at (0.5,0.25) from the pe¬ 
riodicity of the charge density. These small peaks are 
necessary for this structure to exist, and could be used 
to falsify these spin polarons as the origin of the stripes. 


IV. CONCLUSION 

In summary, the microscopic physics of doped cuprates 
has been studied using state of the art quantum Monte 
Carlo techniques. The properties are very close to that 
seen in experiment, with the Fermi surface near that 
of the experiment, and optical excitations in agreement. 
The metal-insulator transition of the cuprates likely oc¬ 
curs between x = 0.00 and x = 0.125, which means 
that this transition has been captured accurately within 
a first-principles method, a new result. In a few years, 
it will likely be possible to use these techniques to study 
the metal-insulator transition in detail. Much of the phe¬ 
nomenology of the cuprates is consistent with the picture 
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FIG. 5: Gharge (a) and spin (b) density in PBEO for two 
holes in a 4x4x1 supercell. The corresponding Fourier trans¬ 
forms for charge (c) and spin (d) have a Bragg peak at 
(0,0.125) because of the periodicity in addition to the plotted 
values. The k values are in units of the reciprocal lattice of 
the primitive tetragonal crystallographic cell. 


emerging from the FN-DMC results. This is a particu¬ 
larly exciting result: since the FN-DMC methodology 
employed here does not use any effective parameteriza¬ 
tion of the interactions, it is a predictive methodology 
that could be helpful in the search for materials with 
similar physics. 

The picture emerging from the FN-DMC calculations 
is that of a spin polaron. The particular polaron found in 
this work was to my knowledge first proposed by Emery 
and Reitei!^ early on, and has been considered by a num¬ 
ber of authors since theiP. The uniqueness of this study 
is that 1) it is truly first principles and explicitly cor¬ 
related, with no adjustable parameters, 2) it was shown 
that the spin polaron can lead to a similar Eermi surface 
structure as seen in experiment, and 3) there are new 
predictions which allow the proposal to be evaluated ex¬ 
perimentally. This spin polaron can be further studied 
using quantum Monte Carlo techniques, and it will be 
fruitful to study the properties of this quasiparticle; for 
example, the interaction between two spin polarons. This 
spin polaron is not stable in the closely related material 
Ca 2 Ni 02 Cl 2 , and thus may be one reason for the unique¬ 
ness of the cuprates. 
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